function result = gauss_int2d_tri(function_f, a, b, c, N_gauss_int2d)

[gauss_points, gauss_weights] = gauss_points_weights_tri(N_gauss_int2d);

det_J = abs((b(1) - a(1)) * (c(2) - a(2)) - (b(2) - a(2)) * (c(1) - a(1)));

sum = 0.0;
for i = 1:N_gauss_int2d
    xi = gauss_points(i,1);
    eta = gauss_points(i,2);
    x = (1-xi-eta)*a + xi*b + eta*c;
    f_val = function_f(x(1),x(2));
    sum = sum + 0.5*gauss_weights(i)*f_val;
end
result = sum*det_J;

end